iT邦幫忙

2024 iThome 鐵人賽

DAY 18
2
自我挑戰組

AI救我系列 第 18

Day 18 - 亥姆霍茲方程-球面波動畫in Python

  • 分享至 

  • xImage
  •  

今天一樣用亥姆霍茲方程去算求面波,其中U(r)波動函數,描述波在空間中的分佈,公式為:

https://ithelp.ithome.com.tw/upload/images/20240930/20168442dOxoLiyCgN.png

接著就來根據先前平面波的程式碼做微調,改成套用求面波的公式:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation


#定義球面波
def spherical_wave(A, k, r, r0, w, t): #A震幅 k波頻 r向量 r0波源的空間座標 w角頻率 t時間
    dr = np.sqrt(np.dot(r-r0, r-r0)) #r跟r0之間的向量差,得出仍為向量。
    #向量和向量做內積平方後取平方根得出絕對值
    phase = complex(0, k*dr-w*t) #k*dr為純量
    #exponential使用負數complex 0為實部,後者為虛部
    return A*np.exp(phase)/dr #根據球面波公式虛除以r

#定義波三維方向波速
k = 2*np.pi

#定義繪圖範圍
D = 5

#繪圖細緻度
N =151

#定義空間座標的2維陣列
x_list = [-D/2 +D*i/N for i in range(N)]
y_list = [-D/2 +D*j/N for j in range(N)]

#np.meshgrid產生網格 製作平面
x_2d, y_2d = np.meshgrid(x_list, y_list)

A = 1 #定義震幅
w = 2*np.pi #定義角頻率

#定義電場為2維的矩陣
E_field = np.zeros((N, N), dtype = complex) 


#定義動畫用的t時間範圍
t_list = [0.1*i for i in range(20)] #跑20個點

#定義畫圖區域
fig, ax = plt.subplots()

#定義r0
r0 = np.array([0, 0]) #只有一個變數,只是該變數是list裡面有兩個數字

#定義動畫更新,將t帶入
def evolution(t):
    ax.clear() #將上一張圖清除
    #每個迴圈都會帶入x和z的數值,return R數值
    for i in range(N):
        for j in range(N):
            r = np.array([x_2d[i,j], y_2d[i,j]])
            E_field[i, j] = spherical_wave(A, k, r, r0, w, t)
            #定義電場 0改為t
    #取電場的實部
    E_real = E_field.real
    #畫圖使用contour (not plot因為plot只能有兩筆資料,三筆資料兩個變數以上就必須用contourf)
    
    ax.contourf(x_2d, y_2d, E_real, cmap = "gray", levels = 99) #cmap為等高線圖色階 levels等高線數 
    
#製作動畫
ani= animation.FuncAnimation(fig, evolution, t_list)
ani.save("Spherical wave_animation.gif")
plt.show()

得出動畫如下:
sphericalwave animation

這次程式的改動有:

  1. 根據球面波公式新增球面波定義,其中新增r0作為球面波源。
  2. 因為是球面波,沒分kx/ky/kz 統一為k。
  3. 定義r0為坐標0,0的波源。
  4. Evolution中改為套用E_field[i, j] = spherical_wave(A, k, r, r0, w, t)。

這樣球面波的動畫就完成啦!


上一篇
Day 17 - 亥姆霍茲方程-平面波(2)動畫繪製in Python
下一篇
Day 19 - 球面波干涉動畫 in Python
系列文
AI救我31
圖片
  直播研討會
圖片
{{ item.channelVendor }} {{ item.webinarstarted }} |
{{ formatDate(item.duration) }}
直播中

尚未有邦友留言

立即登入留言